//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 parenting_age1 parenting_age3 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}


global inputs_std cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std parenting_age1_std parenting_age3_std $outputs

rename cum_avg_daycare_36m_sum c
rename parenting_ages13_std    p
rename parenting_age1_std      p1
rename parenting_age3_std      p3
rename sex male 

global sex0 & male == 0
global sex1 & male == 1
global sex2 & male !=. 

// distribution of inputs and outputs
foreach var of varlist p p1 p3 {
	foreach sex in 0 1 2 {
		foreach sib in 0 1 {
			bootstrap, strata(bwg site tg) reps($bootstraps): reg `var' tg if twin == `sib' ${sex`sex'}, robust 
			matrix b = e(b)
			matrix V = e(V)
			matrix  mc`var'_tw`sib'_sex`sex' = b[1,2]
			matrix  md`var'_tw`sib'_sex`sex' = b[1,1] 
			matrix  se`var'_tw`sib'_sex`sex' = sqrt(V[1,1])
		}
	}
}



// labels
foreach var of varlist p p1 p3 {
	foreach sex in 0 1 2 {
		foreach sib in 0 1 {
			foreach stat in mc md se {
				local `stat'`var'_tw`sib'_sex`sex' = `stat'`var'_tw`sib'_sex`sex'[1,1]
				local `stat'`var'_tw`sib'_sex`sex' = string(``stat'`var'_tw`sib'_sex`sex'',"%9.2f")
			}
		}
	}
}


// parenting
cd $output
foreach var of varlist p p1 p3 {
	foreach sex in 0 1 2 {
		foreach sib in 0 1 {
			# delimit
			global boxp`sib'  text(.375 97.5
				 "{bf: Control Mean}"
				 "{bf: `mc`var'_tw`sib'_sex`sex''}"
				 " "
				 "{bf: {&Delta}}" 
				 "{bf: `md`var'_tw`sib'_sex`sex'' (s.e. = `se`var'_tw`sib'_sex`sex'')}"
			, place(c) box size(medsmall) just(c) margin(l+1 b+1 t+1 r+1) fcolor(none)); 
			# delimit cr

			#delimit
			twoway  (kdensity p if tg == 0 & twin == `sib' ${sex`sex'}, lwidth(dash) lcolor(gs0) lpattern(dash) lwidth(thick)             yline(.45, lpattern(solid) lcolor(gs14)))
					(kdensity p if tg == 1 & twin == `sib' ${sex`sex'}, lcolor(gs6)  lwidth(thick) yline(0, lpattern(solid) lcolor(gs14)) xline(96 , lpattern(solid) lcolor(gs14))),	  
					  

					legend(label(1 "Control") label(2 "Treatment") position(6) size(medsmall) order(1 2) rows(1) region(lcolor(gs0)))
					  xlabel(96[2]102, labsize(medium)        grid   glcolor(gs12) glpattern(solid)) 
					  ylabel(0[.15].45,   labsize(medium) angle(h) glcolor(gs12) glpattern(solid)) 
					  ytitle("Density", size(medium))
					  xtitle("Factor Score")
					  graphregion(color(white)) plotregion(color(white)) ${boxp`sib'};
			#delimit cr
			graph export `var'_tw`sib'_sex`sex'.eps, replace
		}
	}
}
